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Abstract 

We apply the Quasi Monte Carlo (QMC) and recursive numerical integration 
methods to evaluate the Euclidean, discretized time path-integral for the quan¬ 
tum mechanical anharmonic oscillator and a topological quantum mechanical 
rotor model. Eor the anharmonic oscillator both methods outperform stan¬ 
dard Markov Chain Monte Carlo methods and show a significantly improved 
error scaling. For the quantum mechanical rotor we could, however, not find 
a successful way employing QMC. On the other hand, the recursive numerical 
integration method works extremely well for this model and shows an at least 
exponentially fast error scaling. 

Keywords: recursive numerical integration, quasi monte carlo, quantum 
mechanical rotor, anharmonic oscillator, lattice systems, low order couplings 


1. Introduction 

Markov Chain Monte Carlo (MCMC) is the method of choice for simulations 
of quantum field theories or systems in statistical physics. The advantage of 
MCMC is that it can be applied very generally to many physical models. It 
allows to compute expectation values of physical observables (O) with an error 
A which scales only as A oc l/-\/]V, however, where N is the number of samples. 
This error scaling law leads to a very large numerical effort if another significant 
digit in the accuracy of an observable is needed. 

In quantum field theory, in particular quantum chromodynamics (QCD) - 
our theory of the strong interaction between quarks and gluons - very significant 
progress has been achieved in the last years through improvements of the MCMC 
methods used; see, e.g., ref. [1] for an overview. But, even though lattice 
QCD simulations of the theory could be accelerated substantially, computations 
typically run several months or even years on state of the art supercomputers. In 
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addition, none of those improvements have changed the error scaling of 1/ ^fN. 
It will therefore be very demanding to obtain a significant improvement of the 
accuracy of physical observables in this field. 

On the other hand, it is known that Quasi Monte Carlo (QMC) [21 [3] or 
recursive numerical integration mm methods show a much improved error 
scaling. For QMC methods this error scaling reads A oc 1/N°‘ where a can 
reach values of a = 1 or even larger. For recursive numerical integration methods 
the error scaling can be even better and, in some cases, it is even faster than 
exponential, see also the discussion below. 

Clearly, such QMC and recursive numerical integration methods could, thus, 
lead to a much enhanced accuracy of simulations. However, these methods 
have not been tried for generic quantum field theories, so far, and neither their 
applicability nor whether they lead to an improved error scaling is clear. 

In [6l [Zl [8] we initiated a test of QMC methods for the harmonic and anhar- 
monic quantum mechanical oscillator discretized on a Euclidean time lattice; 
see the next section for an introduction to these systems. The results of these 
investigations have been very promising. For the harmonic oscillator, which is 
a Gaussian system, we found an error scaling of A oc 1/A which is optimal. 
Adding a non-Gaussian term in case of the anharmonic oscillator, we found 
a « 0.75 which is not optimal but significantly better than the error scaling of 
MCMC methods. However, it needs to be mentioned that for certain system 
sizes, that is, large Euclidean times, the QMC method did not work particularly 
well for the anharmonic oscillator and no error scaling improvement could be 
established; cf., |7] for details. 

However, QMC methods can also solve another problem: when using MCMC 
for the anharmonic oscillator, only samples in the vicinity of one minimum of 
the action are drawn and the probability to jump to another minimum is very 
small. This leads to a very large, in fact exponential, autocorrelation time. As 
demonstrated in refs. liHi with QMC using the harmonic approximation 
this problem is completely overcome. This is not surprising since QMC avoids 
long autocorrelation times essentially by construction. We stress that, with 
the new developments described in section |5.2[ the iterated numerical integra¬ 
tion method allows to solve the anharmonic oscillator without any problem of 
autocorrelations, as well. 

In this paper, we want to extend the work of refs. oimH] in two directions. 
One direction is to apply recursive numerical integration methods for the an¬ 
harmonic oscillator. As we will see below, this method does not suffer from the 
shortcomings of the QMC method for large Euclidean times. The other direction 
is the investigation of a topological quantum mechanical model, the quantum 
rotor. This model has some characteristic features that are also present for 
non-linear cr-models or gauge theories which are essential and most important 
models for describing elementary particle interactions; see, e.g., refs. laiinKii] 
for introductions to lattice field theories applied to particle theory. However, 
since the quantum rotor is a much simpler 1-dimensional quantum mechanical 
model, it is easier to treat numerically. In this way, it becomes possible to test 
in detail, whether QMC or recursive numerical integration methods work for 
such a model. Clearly, in case the application of one of the methods fails, it will 
become almost impossible to proceed with higher dimensional gauge theories. 
As we will discuss below, so far we have not been able to apply QMC success¬ 
fully to the quantum rotor. On the other hand, recursive numerical integration 
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methods turn out to be highly successful and show an extremely good error 
scaling behavior. 

Another aspect investigated here is the rapid increase of the autocorrelation 
time for the quantum mechanical rotor when the continuum limit is approached. 
For the sake of illustrating this generic problem of MCMC methods, we will 
perform a calculation of the quantum rotor with the Metropolis algorithm. QMC 
and recursive numerical integration methods do not suffer from this problem 
and are, hence, also very advantageous in this respect. Since for the quantum 
mechanical rotor a particular MCMC method which also avoids the problem 
of very large autocorrelation times, the cluster algorithm, can be applied, we 
compare the recursive numerical integration method to this optimal MCMC 
algorithm to see whether a gain can still be found. 

The paper is organized as follows. In the next section we introduce the 
models which are investigated throughout the paper. In section 3 we shortly 
summarize our (failed) attempts to solve the quantum rotor with QMC methods. 
In section 4 we introduce the recursive numerical integration method and discuss 
its theoretical basis. Finally, in section 5 we present our results and conclude in 
section 6. 

2. Lattice systems: 1-dimensional models 

In this paper we investigate two different quantum mechanical models in 
Euclidean tim^ In order to evaluate the models numerically we will discretize 
time and solve the resulting high dimensional integrals numerically. In this sec¬ 
tion we start with some general considerations about calculating observables of 
lattice systems before we describe the two models we investigate, the topological 
oscillator and the anharmonic oscillator. 

Calculating observables in discretized time. The coordinate x{t) € M describes 
the trajectory of a particle in time. Classically, this path of a particle propa¬ 
gating from xo = x(0) to Xi = x(T) during a time period T is determined by 
the minimum of the action S(a;) (with respect to x{t)), 



with the Lagrangian C{x,t) containing all necessary information about the 
model. 

In order to obtain a numerically tractable and positive action while keeping 
the quantum mechanical observables real we have to define the trajectory x{t) 
on a Euclidean, equidistant time lattice with lattice spacing a. Corresponding 
to this, we replace the following quantities by their discretized counterparts: 


t ^ ti := i ■ a, i S {0,1, 2,..., d — 1}, 


( 1 ) 


x{t) Xi := x{ti), 


^Euclidean time is common practice in the standard lattice approach. 
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with T = a ■ d, such that 


T-7 ^ / \ /r,\ 

dt^a} , — ^VXi = -(Xi+i- x^). ( 2 ) 

dt a 

i=0 

We remark that the choice of the discretization of the derivative is not unique 
and alternative discretizations can lead to different error expansions of observ¬ 
ables in terms of the lattice spacing in the continuum limit, a —>■ 0 and T —>■ oo. 
We will use cyclic boundary conditions Xd '■= Xd[ mod d) = throughout this 
paper. 

In a quantum mechanical system not only the classical path contributes to a 
given observable, but all possible paths have to be taken into account. Following 
Feynman’s description, the quantum mechanical system is defined by the path 
integral 


/ e-®Wda;, (3) 

J D<>- 

where Z) is a domain in K. In (§ the transformation to Euclidean time has 
been performed already. For a time discretized quantum mechanical system ([^ 
is well-defined, although it may have a high dimension d, which could be 1000 
or larger. 

The expectation value {0[x\) of an observable 0[x\ = 0 {xq, X 2 , ■ ■ ■, Xd-i) of 
a quantum mechanical model with a discretized action S'[a;] can now be calcu¬ 
lated using the path integral formalism 


(ON) 


JjjaO[x]e 

Jjjd e-^i^Ux 


( 4 ) 


We note that the formalism for treating quantum mechanical systems, as 
it is sketched out above, can be generalized to quantum field theories. Such 
quantum field theories are the actual basis for the theoretical investigation of 
elementary particle interactions. 


Topological oscillator. As indicated previously, the first model we are going 
to consider in this article is the topological oscillator or quantum rotor which 
describes a particle with mass Mq moving on a circle with radius i?, and corre¬ 
spondingly, has a moment of inertia of / = MqR^ . We investigate this particular 
model because it goes beyond the classical quantum mechanical oscillator (de¬ 
scribed later on) and already shows some characteristic features of non-linear 
cr-models and gauge theories which are of prime importance in particle physics. 

The free coordinate of the system is the angle (f) G [—7r,7r), describing the 
position of the particle on a circle with radius i? around the origin. The system 
is described by the action 



and is obtained from the action of a free particle moving in two dimensions, 

S{x)=J + yit)'^) dt (a;,?/)GK^ (5) 
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and the transformation x{t) = Rcos{(j){t)) and y{t) = Rsm{(j){t)) respectively. 
The corresponding discretized action reads 

j d— 1 

= “ X! 

i=0 

where we use the cosine to describe the kinetic part of the action, as ^^(1 — 
cos((/>i+i — 4>i)) = ^ + 0((Vi5ii)^) and hence corresponds to leading order 

to the naive lattice derivative (as in ([^, right). In the following we will, if not 
stated otherwise, set the lattice spacing to a = 1. This means, in particular, 
that the continuum limit is reachecQby T = d • o —> oo. 

One characteristic quantity of the quantum mechanical rotor is the topolog¬ 
ical charge of the system. It describes the number of complete revolutions the 
rotor performs during a time period T 



We use the discretized version 

^ d-l 
i=0 

As an observable, we will investigate the topological susceptibility 

_ (<5^ [(/>]) T^oo I 

^ 47r2J’ 

where {Q‘^[4‘]) is calculated according to @. 

Other important observables of the system are the energy gaps which can 
be extracted from Euclidean correlation functions r(j), 

^d-l 

j (7) 

i=0 

It measures the correlation of angles at different lattice sites separated by a 
distance j. This correlation function has an exponential decay rate with the 
distance j, 

r(j)c5;e-^'^^, j»I. (8) 


From this the energy gap AE between the ground state and the first excited 
state and therefore also the correlation length ^ can be determined. 


c = 


1 

AE 


T-i-o 


21 . 


( 9 ) 


Other energy levels can be computed in a similar way. 


®In principle, the physical extent of time lattice T = d ■ a should be kept constant and 
hence the continuum limit a —>■ 0 requires d —> oo. 
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Through the energy gap, a connection to the topological susceptibility can 
be established, as well; namely. 


Xt = 


T 




AE. 


This, however, only holds in the continuum limit a —>■ 0 and T —>■ oo. 


Anharmonic oscillator. The second quantum mechanical system considered in 
this article is the harmonic oscillator which describes a particle with mass Mq 
moving along a path x G M in a potential proportional to Adding a x^ term 
to the potential, the system is called the anharmonic oscillator. The Lagrangian 
is given by 


£(x, t) 


Mo 

2 



2 


x^ + Ax 


4 


where A G M. In order to keep the action bounded from below, the coupling 
A has to be chosen positive. With A > 0 present, the constant can be chosen 
arbitrarily and for > 0 one finds a distorted harmonic potential while for 

< 0 a double well potential appears. 

Using the discretization scheme mentioned in Q and ([^, the lattice action 
for the anharmonic oscillator becomes 

S[x] = aY^ ^(Vx,)^ + + Axf. 

i=0 

One type of characteristic observables of this system are powers of the position 
X of the particle, i.e., (x), (x^), (x^), etc. 

Another type of observable is again the energy gap which can be extracted 
from the correlator r(j), as already described in eq. 0, and ^ for the topo¬ 
logical oscillator. After exchanging the variables in eq. 0 by the ones of the 
anharmonic oscillator model a similar calculation has to be done here. 

The harmonic and anharmonic quantum mechanical oscillators have been 
studied with QMC methods in iini. 


3. Remarks on Quasi-Monte Carlo attempts 

In this section, we would like to give an overview of our attempts to solve 
the topological oscillator with QMC methods utilizing (randomized) Sobol’ se¬ 
quences (cf., mm)^ The section is intended for researchers who are familiar 
with QMC and may be skipped by non-experts. The addressed attempts were 
motivated by the positive results obtained for the (an)harmonic oscillator model 
(cf., e.g., [5]). At the present state, however, we were not able to observe pos¬ 
itive results with our selected QMC constructions for the topological oscillator 
model. In the following, we will summarize the sampling techniques that were 
tested in combination with (randomized) Sobol’ sequences. 


•^These Sobol’ sequences are low-discrepancy sequences, i.e., they are designed to be more 
uniformly distributed than pseudo-random numbers. 
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Naive Sampling. The naive sampling here simply means that we used the (ran¬ 
domized) Sobol’ points as an approximate uniform distribution and used these 
points directly to evaluate the (re-scaled) integrals according to 



f{x)dx 


1 ^ 


where the Xi are the d-dimensional Sobol’ points. As to be expected, it produced 
good results for small ^ but loses accuracy and convergence speed as ^ increases 
since most samples have little to no weight rendering them irrelevant. In fact, 
as ^ grows significantly larger than | we have not been able to observe any 
convergence in the range of sample sizes that are feasibly generatable. 


Harmonic Sampling. Since we observed good results using the Sobol’ sequences 
with the harmonic and anharmonic oscillator, it seemed reasonable to approxi¬ 
mate 


cos {(j)i - (j)j) «!-(())*- (jijf ', 


i.e., to draw just as in and re-weight to account for the error made in drawing 
from this approximate distribution. Using the cumulative (standard) normal 
distribution <i>, there are three obvious approaches to inverting $ with range 

[-7r,7r). 

• Using K as a covering space of [—tt, tt) yields the “inverse” 


where 
MOD : 


(0,1) —>■ [—TT, tt); a; ^(x) MOD 27r 


X 1 
^ 2 


[—TT, tt); X ^ X — round ^ —^ 27r = a: — 


27r. 


• Use the complete normal distribution and dismiss points not in [—7r,7r]‘’*. 
(Note that this approach destroys the property of uniformity of the Sobol’ 
sequence by rejecting some of the points.) 

• Restrict to [— tt, 7r)-slice of the normal distribution, i.e., invert 


$(a;) 


$(a;) — <i)(— tt) 
2$(7r) - 1 


Qualitatively, all three approaches yield the same result; our quadrature was 
highly instable yielding seemingly random numbers and we have not been able 
to stabilize them. We think this is due to an inherent under-representation 
of samples in the region where \4>i+i — 4'i\ is close to 27r. These samples have 
a large weight but the chance of drawing them is slim. Hence, hitting these 
regions a little more often can make a significant difference which we observe as 
an instability of our results. 
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Inversive Samplings. Inversive samplings (cf., [H] chapter 11.2) go one step 
further than the harmonic sampling by choosing better approximations than 
the harmonic one. However, we still have the constraint that we need to be able 
to effectively draw from the distribution. In this case, we are looking for an 
expression of the form 

d-l 

exp(-S'[(/)]) « Y[pj{<l)j)- 
j=o 

In fact, the two main examples we tested were of the form 

d-l 

exp(-S'[(/)]) « YIp{(I)j) 
j=o 

where p is a polynomial or a step-function. The step-function has the advan¬ 
tage that we can draw with very little effort. However, QMC does not like 
discontinuous integrands, i.e., we expect to lose convergence speed. Choosing 
a polynomial p does not have that draw-back but generating samples is a lot 
harder since it involves numerically inverting the cumulative distribution func¬ 
tion defined by the polynomial used as a density. 

Our results are satisfactory in the sense that we did not observe an increase 
in autocorrelation time going to small lattice spacings (in fact, we would be 
highly surprised if that happened since we are not using a random walk to 
generate samples) but we have not been able to improve over standard Monte 
Carlo methods where they are applicable, that is, to obtain an error scaling 
better than O where N is the sample size. 

Sampling - Importance Resampling. SIR (cf., [HI [161 [17] ) is a closely related 
concept where a pool of sample points is generated (here, using (randomized) 
Sobol’ sequences) and stored. Samples are then to be drawn from the discrete 
distribution of points in the pool with respect to their weight in the target 
distribution (cf., e.g., [Hj chapter HI.2). It worked fairly well for small sample 
sizes but, trying to go to large numbers of samples, the pool must increase as 
well in order to keep the systematic error of the changed distribution small 
rendering this method too memory demanding to be viable for high accuracy 
calculations (an error analysis of QMC based SIR methods can be found in [18jl. 

Envelope Inversive Rejective Samplings. An approach to draw directly from the 
target distribution but with the help of inversive sampling starts by choosing 
an approximation S');/)] of such that 

e-SW < 

We call such an inversive sampling an Envelope Inversive Sampling. Drawing 
from the envelope, we may now use a rejection step (cf., [H]) to ensure drawing 
from the target distribution in the usual Monte Carlo manner. 

Our results are comparable to the Inversive sampling results, that is, no 
increased auto-correlation for small lattice spacings but also no accelerated con¬ 
vergence with respect to standard Monte-Carlo methods. 


Smoothed Envelope Inversive Rejective Samplings. Here, the idea is that we 
might obtain an edge starting from the envelope inversive rejective sampling 
by smoothing the integrand; viz., the rejective step is a step-function in the 
integrand which behaves poorly with QMC and, hence, may yield better results 
if subjected to a smoothing operation. The method of Smoothed Rejection can 
be found in m and effectively replaces the step-function describing the rejection 
step by a continuous function and additional re-weighting. However, we could 
not observe any improvement over the envelope inversive rejective sampling. 


4. The method of recursive numerical integration 


In this section, we consider the method of recursive numerical integration 
(see HH] and references therein) also sometimes called the method of iterated 
numerical integration, for the approximation of the integration problems in 1- 
dimensional lattice theory stated in section The main idea of this method 
is based on the fact that if the integrand at hand can be described as the 
product of low-dimensional functions, each one describing the interactions or 
couplings between few consecutive objects, then we can write the final high¬ 
dimensional integration problem as the iteration of coupled low-dimensional 
integrals. At this point a quadrature rule for low-dimensional integration has to 
be chosen in order to carry out the successive approximation of the underlying 
low-dimensional integrals. In the following, we will consider the case of simple 
1-neighboring couplings for simplicity. For a description of the general case with 
several neighboring couplings or branchings, we refer again to [S]. Thus, we have 
an integrand function of the form 


d-l 


fix) = n Mxi,X(^i+i)(d)) 


2 = 0 


and we would like to approximate 

^ d-l 


I = 


! D'^ 


;=o 


where the notation (.)(d) means to take the argument modulo d, and H C K is 
a 1-dimensional domain. Thus, by the Fubini-Tonelli theorem, assuming that 
the desired integral value exists, we can write the problem equivalently as 

^ ^Id ifjydoi^o,x:i)... 

■' U fd-3{^d-3,^d-2) (^J fd-2{^d-2,Xd-l)fd-l{oi:d-l,^o)dXd-l^ dXd- 2 '^ ■ ■ ■ 

... dxi^dxQ. 

Furthermore, we can consider scalings co,...,Cd-i > 0, and define I* := 
(nr=ro^ ir) d, such that 

- f f f foixo,xi) 

JD ^Jd Co 

fd-3{xd-3,Xd-2) f f fd-2(xd-2,Xd-l)fd-l(xd-l,XQ) \ ^ \ 

ID Cd-3 \J D 


r = 


Cd-3 


Cd-2Cd-l 


... dxi^dxf). 
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Usually the quantities cq, , Cd-i will have to be chosen adaptively as they are 
used in order to avoid under/over-flows in the recursive method for high dimen¬ 
sions due to limited machine accuracy. By selecting an adequate 1-dimensional 
quadrature rule Q with m points and weights for approximating the underlying 
integration problems over D in each iteration, we can write the iteration method 
in recipe form as 

1. Fix, if possible, a quadrature rule with m points x^,... and weights 
wi,..., Wm that works well for one dimensional integrals of type 


[ fi-i{0i,z)fi{z,92)dz, for 1 < t < d- 1 , 6 »i, 6»2 e ( 10 ) 

Jd 


2. Use the quadrature in 1. to estimate 

1 


Fd_i{xd-2,Xo) ■■= - ^ - / fd-2{Xd-2,Xd-l)fd-l{Xd-l,Xo)dXd-l 

Cd-2Cd-i Jd 

^ m 

- y^,Wjfd-2ixd-2,X^)fd-l{x^,Xo), 


Cd-2Cd-l 


i=i 


over the grid of points {x ^,..., x™} x {x^,..., x™} C By defining the 
transfer m x m matrix 

•.= fi{x^,x^), l<k,l<m, 0<i<d—l, 

we can write in matrix form 

1 


[Fd-l{xd-2J xy)\l<ij<r 


Cd-2Cd-l 


-Md- 2 diag{{wi,.. .,Wm))Md-i, 


where 


diag{{wi,.. .,Wm)) ■= 


(wi 

0 • • 

■ M 

0 

W2 ■ ■ 

0 


0 • • 

• Wra) 


3. For i = d — 2,d — 3,... ,2 approximate iteratively 

Ui(xi_i,xo) = ^— / /i_i(xi_i,Xi)Fi+i(xi,xo)da 

Ci-i Jd 
1 


Q— 1 


Wj-/i_i(Xi_i,X-^)fi+l(x^Xo) 




over the grid of points {x^,... ,x™} x {x^,..., x'"} C D"^. Thus, the result 
of step 3. over {x^,..., x™} x {x^,..., x™} can be written as the matrix 
products 


[F 2 


d-2 ^ ^ 

( II —Midiag{{wi,.. .,Wm))) - Md-i, 

' • 1 Cf / Cd—l 

t—1 


with the notation for matrix products ■.= B^B2...Bd-2. 
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4. Estimate with the m-point quadrature in 1. the function 


n(x„,x.)= [ 

Jd Co 

over the grid {xq, ... ,x^} C D. Finally, estimate with the m-point 
quadrature in 1. the function 


I* = Fi{xo,xo)dxo. 


'D 


The result of this step can be written in matrix operations as 


/ d-2 ^ 

* « Tr I diag{{wi,.. ■, Wm))(| —Midiag{{wi,.. .,Wm))^ - Md-i 

\ i=0 ^d-l 


d-1 


= Tr —Midiag{{wi ,..., w^)) . 


\i=0 


As mentioned in the previous section we are mainly interested in consider¬ 
ing two types of integrands for high-dimensional integration. The first type of 
integrand is given as the weight of the Boltzmann distribution. In this case we 
usually obtain that the functions /i,0<i<d—1, satisfy /o = /i = • • • = fd-i 
due to isotropic conditions of the model. Thus, this case yields Mq = Mi = 
■ ■ ■ = Md-i =: M and we can choose Cq = Ci = • • • = Cd-i ='■ c. Hence, the 
identity above reduces to 


/* = Tr 



Mdiag{(wi, 





Note that the similarity relation 


Mdiag{{wi,Wm)) ~ diag{{y/wi ,..., y/wm))Mdiag{{y/w {,..., 

holds, where for m x m matrices A, B we say that A ^ H if and only if 
A = CBC~^, for an invertible matrix C. The resulting matrix on the right 
hand side of the relation above is symmetric and efficient algorithms are known 
for its diagonalization (cf., [H]). By use of the eigenvalue decomposition of a 
diagonalizable matrix A, we can write = CD‘^C~^ where D is the diagonal 
matrix of eigenvalues, and C is the corresponding invertible matrix of eigen¬ 
vectors. Because the trace of a diagonalizable matrix equals the sum of its 
eigenvalues, to calculate the trace of A^* we only need to calculate the eigen¬ 
values of A, raise them to the power of d, and finally sum them up. It is also 
worth to mention that for particular applications, recursive-multiplication for 

calculation of {^Mdiag{{wi ,..., Wm))^ may be also a competitive method. 
The second type of integrand is given by the Boltzmann weight times an ob¬ 
servable function. For simple observable functions usually considered in lattice 
theory, similar gains over direct matrix-matrix multiplications can be obtained 
by condensing a sequence of matrix-matrix multiplications into a power form. 
Particular examples and the corresponding implementations will be discussed 
in detail in the following section. 
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5. Numerical experiments 

In this section, we will investigate the applicability of the method of recursive 
numerical integration to the topological rotor and the anharmonic oscillator, 
as described in section In particular, we will compare the calculations of 
recursive numerical integration methods with standard MCMC methods applied 
to lattice systems dn Eg and also with improved methods, i.e., the cluster 
algorithm [53] and randomized QMC methods for the anharmonic oscillator [7] . 

5.1. Topological oscillator 

For the sake of illustrating the difficulties that can appear in standard 
MCMC calculations, we will perform initial computations with the Metropolis 
algorithm [24|. This algorithm, although very generally applicable, is often not 
the optimal choice to simulate a model and is usually replaced by a more suit¬ 
able technique. However, here we use the Metropolis algorithm to demonstrate 
the generic difficulty of many MCMC methods that towards the continuum 
limit the autocorrelation time grows rapidly. At fixed number of samples, this 
leads to very large errors when the lattice spacing is reduced and, correspond¬ 
ingly, to highly increasing computational costs of the simulations. Throughout 
the discussion below, we use a moment of inertia of / = 0.25 for the quantum 
mechanical topological rotor. 

Besides promising a much improved error scaling, recursive numerical inte¬ 
gration methods avoid this explosion of the autocorrelation time and are there¬ 
fore highly superior to MCMC techniques. Nevertheless, for special situations 
there exist MCMC methods which also avoid the increase in autocorrelation 
time. For the topological rotor investigated here, cluster algorithms [53] can be 
applied which exhibit a basically constant behavior of the autocorrelation time 
with respect to the lattice spacing while still showing the l/-\/iV error scaling 
behavior of MCMC methods. Hence, we will explore how the recursive numer¬ 
ical integration technique compares to an optimal MCMC method such as the 
cluster algorithm and whether a gain can still be found. Since for the quantum 
mechanical topological rotor we could not find a satisfactory realization of QMC 
methods to solve the system, we will not discuss this approach in this section. 

5.1.1. MCMC method 

A basic algorithm of the MCMC class is the Metropolis algorithm [24] which, 
as discussed above, will be used for illustration purposes only. The Metropolis 
algorithm uses importance sampling as described in section]^ to create samples 
of a given system on which observables O can be calculated. These observ¬ 
ables are then averaged over all samples to give the expectation value (O) with 
an error A(0) oc where N is the number of samples. The problem of 

the Metropolis algorithm - and MCMC methods in general - is that the sam¬ 
ples are not all independent of each other but correlated. This correlation is 
measured through the autocorrelation time which is a property of the employed 
algorithm and, if the so-called integrated autocorrelation time is taken, also of 
the considered observable. 

The behavior of the error and the integrated autocorrelation time of the 
topological susceptibility when applying the Metropolis algorithm as a function 
of the lattice spacing a is shown in figure We have chosen the topological 
susceptibility as an observable, since it is very sensitive to correlations between 
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Figure 1: Logarithmic error and logarithmic integrated autocorrelation time 
behavior of the Metropolis algorithm for fixed number of samples {N = 10®) 
dependent on the lattice spacing a. Left: Error of the topological susceptibility 
Aytop, right: the corresponding integrated autocorrelation time Tint- 


samples and often shows the largest integrated autocorrelation time in a simu¬ 
lation. It is also for this reason that the quantum mechanical topological rotor, 
where topological effects can be studied, is a good test case to compare MCMC 
and QMC or recursive numerical integration methods. 

We remark that for all values of the lattice spacing shown in fig. a fixed 
number of samples (N = 10®) has been used. In addition, for the error cal¬ 
culation the integrated autocorrelation time has been fully taken into account. 
As can be observed in fig. the integrated autocorrelation time Tint (right 
plot) and therefore the error of the topological susceptibility Aytop (left plot) 
are rapidly growing with decreasing lattice spacing a. This leads to highly 
increased computational costs. 

For systems in statistical physics or quantum field theory the continuum 
limit is reached by approaching a critical point. Since the rapid increase of 
the autocorrelation time in this limit is generic for many MCMC algorithms, 
this constitutes a most severe problem when higher dimensional systems are 
explored. 

For the explicit model considered here, there exists a cluster algorithm [23] 
which avoids the problem of an increasing autocorrelation time for shrinking 
lattice spacings. Figure]^ shows the behavior of the autocorrelation time (right 
plot) and the error of the topological susceptibility (left plot) for the cluster 
algorithm as a function of a. As can be seen for the cluster algorithm, there is 
even a tendency that both, the error and the autocorrelation time, shrink for 
smaller a. 

To compare the two algorithms directly, figure shows the observable Xtop 
dependent on a for the Metropolis and the cluster algorithms. For larger values 
of a the error from both algorithms are similar. However, for small values 
of a the error of the Metropolis algorithm becomes so large that it would be 
very demanding to reach an accurate result below a certain value of the lattice 
spacing, thus making a continuum extrapolation of Xtop rather difficult. On the 
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Figure 2: Error and integrated autocorrelation time behavior of the cluster al¬ 
gorithm for fixed number of samples {N = 10®) dependent on the lattice spacing 
a. Left: Error of the topological susceptibility Aytopt right: the corresponding 
integrated autocorrelation time Tint- 


other hand, for the cluster algorithm the error stays practically constant; hence, 
allowing us to reach very small values of the lattice spacing giving a much better 
control of the continuum limit. 

We stress again that the discussion above is only intended for an illustration 
of the generic behavior of MCMC methods and to demonstrate the problem of 
an increasing autocorrelation time when approaching the continuum limit for 
certain classes of MCMC algorithms. As the examples show, MCMC methods 
often have to face the difficulty of very large autocorrelation times which would 
be absent for QMC or recursive numerical integration methods. In the model 
considered here an optimal algorithm can be used, the cluster algorithm. It 
will therefore be interesting to see, whether the recursive numerical integration 
method is still advantageous even for cases when highly improved MCMC tech¬ 
niques can be employed. Although not relevant for this paper, we would like to 
mention that cluster algorithms are not applicable to gauge theories, so far. 

5.1.2. Recursive Gaussian quadrature 

In order to demonstrate the precision that can be reached with the recursive 
numerical integration method even at very small values of the lattice spac¬ 
ing, we will first discuss some numerical results that we have obtained with 
this approach. To this end, we have implemented the method of recursive nu¬ 
merical integration described in section using Gauss-Legendre mesh points 
(abscissae), and applied it to the topological oscillator. Fixing the number of 
integration points to to = 120 and the value of the moment of inertia I = 0.25, 
we calculated the topological susceptibility Xtop, the energy gap AE and the 
ratio of both observables as a function of the lattice spacing; cf., fig. 

For the employed value of / = 0.25, there are theoretical predictions for the 
observables we consider here, see section In particular, in the continuum 
limit we should find for the topological susceptibility Xtop —> , for the 

energy gap AE —^ = 2, and for the ratio —>■ 
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Figure 3: Continuum extrapolation of the topological susceptibility Xtop calcu¬ 
lated by the Metropolis and the cluster algorithms for a fixed number of samples 
{N = 10^). 


In figure we show our results obtained with recursive Gauss-Legendre 
quadrature with the expected continuum values, as given above, subtracted. 
The graphs nicely show that for all observables the results converge to zero; 
thus, being fully consistent with the theoretical expectations. Note that in the 
left column of fig. we use a linear scale for the observables considered while 
in the right column a logarithmic scale is used which demonstrates the high 
precision we can reach with recursive Gauss-Legendre quadrature. Furthermore, 
note that the ratio has a peak at higher values of a, due to peculiar lattice 
artifacts, which is consistent with previous studies of this model [5S]. 

Let us now turn to the interesting question of the error scaling for the re¬ 
cursive integration method. The behavior of the error for the topological sus¬ 
ceptibility is shown in figure for a fixed lattice spacing a = 0.4 as a function 
of the number of integration points m. Note that we use a logarithmic scale for 
plotting the error of the topological susceptibility. We define the error by the 
difference of Xtop obtained for m = 560 and ytop computed at the given number 
of integration points m < 480. For 220 < m < 480 we fit an exponential function 
of the error in m which appears as a straight line in fig. where a logarithmic 
scale is used. The good agreement between the data and this exponential fit 
suggests that asymptotically the error scales down at least exponentially fast. 
Recall that the error behavior for Gauss-Legendre quadrature with m points for 
infinitely differentiable functions scales like ~ where 

the latter relation holds due to Stirling’s approximation formula. 

It is clear that with asymptotic exponential (or even better) error scaling 
the numerical recursive integration method will outperform any algorithm that 
shows an algebraic error scaling, in particular the I/Vn behavior of MGMG 
algorithms. Still, it is an interesting question whether for small, more practical 
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Figure 4: Gaussian quadrature - continuum extrapolation of three different 
observables, the topological susceptibility Xtop (top row), the energy gap AE = 
El — Eq (middle row) and the ratio (bottom row) with m = 120 integration 
points. Each observable is presented with a linear axis scale (left column) and a 
logarithmic axis scale (right column). For all results shown, we have subtracted 
the theoretical value of the given observable in the continuum such that they 
should converge to zero in the continuum limit. 


values of JV (or to) a gain can be obtained from the recursive integration method. 

We, therefore, measured the run-time of the Gauss and the cluster algo¬ 
rithms needed to obtain a given error of the topological susceptibility. For these 
measurements we ran both algorithms with the same input parameters (we used 
here a lattice spacing of a = 0.1) on the same stand-alone computer. We then 
varied the number of samples and integration points for the cluster and recur- 
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Figure 5: Gaussian quadrature - error behavior of Xtop in dependence of the 
number of points m used in the integration and with fixed lattice constant 
a = 0.4. Note that we use a logarithmic scale to exhibit the error. The blue line 
is an exponential fit of the error for 220 < m < 480 which appears - through the 
use of a logarithmic scale - as a straight line for a function of m. The consistency 
of the data with the linear behavior suggest that the error scales down at least 
exponentially fast. 


sive Gauss-Legendre quadrature, respectively, and measured the run-time. We 
repeated each measurement ten times and averaged over the measured times to 
get an error estimate of the run-time. This error originates from the number 
and kind of other processes that are running on the machine at a given time 
and is noticeable for small run-times. For the cluster algorithm, in addition, 
the size and distributions of the generated clusters can vary leading to different 
run-times of the algorithm. 

Being an MGMG method, for the cluster algorithm the procedure of repeat¬ 
ing all runs ten times allows us to also determine the error on A^top, he., the 
error of the error. The recursive Gauss-Legendre quadrature, on the other hand, 
is purely deterministic and, therefore, gives the same result for Xtop without any 
error every time we run it. Since for this test we have used different parameters, 
we have also chosen a different gauge value of m = 400 and calculate the error 
for Xi(m < 400) via the difference to this value, Axi = \xi — Xi{^ = 400) |. The 
error of this error is neglected here. 

In figure the run-time t needed to achieve a given error of the topological 
susceptibility is shown for both algorithms. Note that the graph is shown in a 
double logarithmic scale. The cluster algorithm shows the for MGMG methods 
typical l/y/i behavior, indicated by the red line. The recursive Gauss-Legendre 
quadrature has a much steeper negative slope, such that the calculation to 
reduce the error by a specific amount is substantially faster than using the cluster 
algorithm. Additionally, even for small run-times the recursive Gauss-Legendre 
quadrature shows a much reduced error already. Therefore, the recursive Gauss- 
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Figure 6: Comparison of Cluster algorithm and recursive Gauss-Legendre 
quadrature. We show the run-time t in seconds needed to get a given error 
on the topological susceptibility with fixed lattice spacing a = 0.1 on a stand¬ 
alone computer. Note that the graph is plotted in a double logarithmic scale. 
We used N = 10^... 10® samples for the cluster algorithm and m = 10...300 in¬ 
tegration points for the Gauss algorithm. To achieve the error of the error the 
numerical experiment has been repeated ten times for a fixed set of parameters. 
The run-time error of both algorithms results from different execution times of 
the program due to state of the computer when the program was executed. For 
the cluster algorithm the error in the run-time is also determined from differ¬ 
ent distributions of clusters generated. The red line shows the expected 1/Vi 
behavior of the cluster algorithm as a MCMC method. 


Legendre quadrature is not only superior in the asymptotic regime at large run¬ 
times where it shows a much improved error scaling but is also advantageous 
when only small run-times are employed. This makes it a promising approach 
if one thinks of systems in higher dimensions where only short run-times can be 
afforded. 

5.2. Anharmonic oscillator 

Applications of QMC methods to the anharmonic oscillator model and com¬ 
parison with MCMC techniques have been studied in IS1I711H]. In particular, 
the superiority of randomized QMC techniques applied to the anharmonic os¬ 
cillator model has been observed for small time periods T < 1.5 independently 
of the size of lattice spacing a. Here we investigate whether the recursive in¬ 
tegration method can also be successfully applied to this model and whether 
improvements over randomized QMC can be observed. To this end, we carry 
out tests by choosing different pairs of problem dimension and spacing (d, a). 
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The time period T associated to a pair (d, a) is given by T = da. The choice 
of the pairs (d, a) has been done with the aim to show that the anharmonic 
oscillator model can be approximated satisfactorily with a particular recursive 
numerical integration method. The observables investigated are (x^) and (x'^). 
The rest of the parameters of the model have been fixed for the experiments 
to Mq = 0.5, = —16.0, and A = 1.0. We perform recursive one-dimensional 

integration based on Gauss-Legendre quadrature. To this end, we first re-scale 
the integration variables and re-arrange the terms in the action such that (using 
the same notation as in section]^ we end up with a transition function 


fi ) 






The quartic negative term dominates the quadratic positive term 

— ^^^xf when \xi\ > for i = l,...,d. Outside of the region [ — 


4Aa ■ 


Molii!!' 
4Aa . 


, the term inside the exponential in the transition function 
fi{xi, Xi+i) remains always negative, and decays mainly with a quartic rate that 
adds to the quadratic coupling —{xi — Xi+i)^. Therefore the marginal tails of 
this two-dimensional transition function decay faster than the marginal tails of 
a bivariate normal, and this fact can be used to select a region of main impor¬ 
tance for the one-dimensional parametric integration problems described in (101. 
Thus, we additionally search to ensure that the decay obtained by the dominant 
quartic term overshadows the contribution of the remaining coupling quadratic 
term —(xi — Xj+i)^ in the action. As an heuristic, we then chose the main 

importance region for integration to be of the form [ — ) 

for a number p > 1, and validate numerically that the remaining integration re¬ 
gion ( — oo, U \p\l, -l-oo) has a negligible contribution to the 

integration problem. Note that this is the same as to say that the original inte¬ 
gration problem over can be truncated satisfactorily for recursive numerical 

integration to the integration region [—p\f^i)^,P\f ^ 4 x^^ ]^- The selection of 
p > 1 can be taken to be a positive value that ensures — ^j^xj < —K, 

for some positive value K, and Xi in ( — oo, —p\J ^fxa U \p\J ^fx^ , -l-oo). 
By integrating the tails of a normal density outside of the region 


|t| < 7(7, we obtain a value below (but close to) 10 The maximal value of 


e 

-10 


the exponential function in the density is given by at t = la. Thus, for our 
tests we may select the conservative value of AT = 24.5 based on the fact that 
the tails of the one-dimensional integration problems in the recursive numerical 
integration method decay faster than a normal density, and that inside of the 
main importance region the behavior of the marginals seems mainly similar to 
the one of a shifted normal density with o' = ^ (due to the coupling quadratic 
term). This choice also makes sense since our computations are carried out in 
double precision and we aim to obtain at most 10 digits accuracy for the inte¬ 
gration problems. Thus, at the end, we select p > 1 to be the positive number 
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= -24.5. 
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We choose to take sample sizes m inside of the selected importance region to 
be small multiples of , since we would like to have good integra¬ 

tion accuracy in each unit-length interval inside of the parametric problem (10) 


in the region [ — ■ Note that in terms of the spacing pa¬ 

rameter a, we have m ~ 0{a~i). The remaining marginal tails outside of 

[ ~ estimated with very few Gaussian points (or 

no points at all) if the relative contribution of the tail of the integrand to the 
total integral is negligible relative to our maximal target accuracy. For our ex¬ 
periments, we used Gauss-Hermite points (corresponding to a weight function 
of the form , c > 0 ) and Gaussian points generated from a weight function 

4Aa^ 4 / / 

e “o ‘, for the integration region ( — oo, —p\J ^^ ] U [py , -l-oo). By 

increasing the sample sizes in this region we were able to observe that the rela¬ 
tive contribution of the integrand on \ [ — p\l^^^x^iP\l ^la ^ total 

integrals on seems less than 10 “^^. Therefore we believe that the selected 
importance region for integration can approximate the original problem with 
great relative accuracy (close to 10 “^^). 

Figur[^shows our calculations of the two observables and (cc^), depend¬ 
ing on the sample size m, for different values of lattice points and spacing (d, a). 
Note that the highest m « 4600 was taken for (d, a) = ( 2 ^^, 2 “^^), but even 
in this case the computations took no longer than 3 minutes with a standard 
PG. For high values of T (like T = 4096 ) we observe very good results, as well, 
contrary to what has been observed with randomized QMG in [5]. We would 
like to remark that for the ground state energy Eq , which, by virtue of the virial 
theorem, is related to (a;^) and (x^) by Eq = the result¬ 

ing estimates for T = 4096 matches with the theoretical value in 5 significant 
digits, Eq = 3 . 8636669 , calculated in [ 26 ], namely Eq = 3.86367053759882 for 
(d, a) = ( 2 ^ 4 , 2 - 42 ) ^ _ 4gQQ^ 
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Figure 7: Shown are the estimated values of {x^) (left) and (right) in the 
anharmonic oscillator model obtained with recursive Gauss-Legendre quadra¬ 
ture for the pairs (d, a) = (2’^, 2“’^), (2^^, 2“^), (2^^,2“^°) and (2^"^, 2“^^), with 
corresponding time periods T = 1.0, T = 16, T = 16 and T = 4096. In all 
cases the relative accuracy achieved with the corresponding maximal samples 
m seems to be more than 10 significant digits. 
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6. Conclusions 


In this paper we have applied the methods of Quasi Monte Carlo and recur¬ 
sive numerical integration to two quantum mechanical models discretized on a 
Euclidean time lattice within the path integral approach. 

The first model we have considered is a quantum mechanical rotor which 
is - to some extent - similar to higher dimensional spin systems such as non¬ 
linear tr-models. For the quantum mechanical rotor we could, unfortunately, 
not find a successful implementation of the QMC method to solve this model. 
On the other hand, the method of recursive numerical integration led to a much 
improved accuracy even when compared to an optimal MCMC algorithm for 
which we have chosen the cluster algorithm. Conceptually, the error scaling of 
the recursive numerical integration is at least exponentially fast in the number of 
integration points m. In fig. we could indeed demonstrate that asymptotically 
this exponential error scaling is realized. Figure]^ shows that even for a small 
number of integration points the accuracy of the recursive integration method is 
already much higher than the one of the cluster algorithm with a correspondingly 
small number of samples. 

We remark that we have chosen the cluster algorithm as the MCMC method 
since it avoids the increase of the autocorrelation time towards the continuum 
limit, a feature which is shared by the recursive numerical integration technique 
by construction. The avoidance of autocorrelation times and our finding that 
a very high accuracy can be reached already for a small number of integration 
points makes the recursive numerical integration technique a very promising 
method for more difficult systems where only a small number of samples can be 
realized, e.g. in higher dimensions. 

The other model we looked at is the anharmonic quantum mechanical os¬ 
cillator. Here, we had found earlier that the QMC method shows an improved 
error scaling HEHH]. When the extent of the time lattice is kept short, T < 1.5, 
both, QMC and recursive numerical integration show a comparable performance 
with an improved error scaling which is faster than l/y/N. However, as noted 
in [3 QMC becomes inefficient when the time extent T is made larger than 
T = 1.5, independent of the value of the lattice spacing a. 

Therefore, we tested the performance of the recursive numerical integration 
method for various choices of the dimension d and lattice spacing a. As fig. 
shows, we can choose a very broad range for [d, a) where we still find an ex¬ 
tremely good performance of the recursive numerical integration method. In 
fact, the values of the time extent, T = da can assume very large values such as 
T = 4096 and the values of a can become tiny, e.g., a = 2“^^, while still a rapid 
convergence of the considered quantities is observed. 

In conclusion, we have tested two methods to evaluate the path-integral in 
Euclidean time for quantum mechanical systems. These are the QMC and the 
recursive numerical integration methods. While we could not find a successful 
implementation for QMC in the case of the quantum mechanical rotor, the 
technique of recursive numerical integration has been highly successful with an 
exponentially fast error scaling. It will be very interesting to test this method 
for more complicated 1-dimensional models and, of course, for systems in higher 
dimensions. 
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